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The application of newly developed first-principle modeling techniques to liquid water deepens our 
understanding of the microscopic origins of its unusual macroscopic properties and behaviour. Here, 
we review two novel ab initio computational methods: second-generation Car-Parrinello molecular 
dynamics and decomposition analysis based on absolutely localized molecular orbitals. We show 
that these two methods in combination not only enable ab initio molecular dynamics simulations 
on previously inaccessible time and length scales, but also provide unprecedented insights into the 
nature of hydrogen bonding between water molecules. We discuss recent applications of these 
methods to water clusters and bulk water. 
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I. INTRODUCTION 

Liquid water is of paramount importance for life on 
Earth. That is why its properties and behaviour has been 
a subject of scientific investigation for several centuricsA. 
It has long been established that liquid water has numer- 
ous unusual properties and exhibits anomalous behaviour 
for a wide range of different conditions. Nevertheless, de- 
spite extensive research, the underlying physical origins 
of many of these phenomena remain unclear. 

An isolated water molecule in the gas phase obeys a 
simple geometry and electronic structure. However, in 
the liquid phase, each water molecule forms multiple re- 
markably strong hydrogen bonds (HBs) with its neigh- 
bours and, thus, becomes a structural unit of the ex- 
tended HB network. The complex structure, energet- 



ics and dynamics of this fluctuating network determine 
many anomalous macroscopic properties of liquid water. 
Therefore, a detailed investigation of elementary molec- 
ular processes in the HB network - vibrations, reorien- 
tations, diffusion, HB-rcarrangemcnts - is crucial for un- 
ravelling water mysteries, for better understanding of its 
role in nature and, consequently, for its better utilisation 
in artificial applications. 

Experimental investigations of fundamental processes 
in the HB network at the molecular level has become 
possible only in the last few decades with the advent 
of sophisticated spectroscopic probes with femtosecond 
time resolution. However, even the most advanced spec- 
troscopic techniques measure only the spectroscopic re- 
sponse of the complex system of interconnected water 
molecules, which is often difficult to interpret in terms of 
the real-time molecular structures and rearrangements. 
A recent highly controversial interpretation of the X-ray 
spectra of liquid water as evidence for its "chains and 
rings" structure^ is perhaps the most dramatic illustra- 
tion of the intrinsic ambiguities of spectroscopic analysis. 

The inconclusivcness of spectroscopic measurements 
has made computational modelling an indispensable 
tool for obtaining intelligible information from intri- 
cate experimental data. Ab initio molecular dynamics 
(AIMD)2r— , in which the interactions are obtained from 
accurate electronic structure calculations, has become 
widespread in computational studies of liquid water— ~—. 
However, despite constant advances in high-performance 
computing, the great computational cost of AIMD still 
imposes severe constraints on the length and time scales 
attainable in AIMD simulations. Because of such restric- 
tions the results of calculations often contain finite-size 
errors and/or statistical uncertainties due to insufficient 
sampling, the magnitude of which is difficult to estimate 
accurately. Here, we review a novel computational ap- 
proach to AIMD that accelerates simulations with hun- 
dreds of water molecules such that trajectories as long 
as a nanosecond can be generated. We demonstrate that 
this new approach, which combines the advantages of 
both Car-Parrinello and Born-Oppenheimer molecular 
dynamics (MD), allows us to predict the properties of 
liquid water which are difficult to obtain with less effi- 
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cient conventional AIMD methods. 

The thermodynamics and kinetics of elementary pro- 
cesses in the HB network of liquid water are inextricably 
connected to the strength of interactions between water 
molecules. Therefore, in addition to computationally ef- 
ficient algorithms, there is considerable interest in devel- 
oping accurate electronic structure methods, which pro- 
vide physical insight into the nature of interactions that 
determine the strength of hydrogen bonding. In this re- 
view, we discuss a recently developed energy decompo- 
sition technique, based on absolutely localized molecular 
orbitals (ALMOs), and its applications to water clusters 
that reveal an interesting and somewhat unexpected view 
of the electronic origins of hydrogen bonding. We also 
show how the combined application of the novel AIMD 
and energy decomposition methods help to address con- 
troversial questions related to the symmetry and order of 
the HB in liquid water. 



II. COMPUTATIONAL METHODS 

A. Second-generation Car-Parrinello molecular 
dynamics 

Until recently, ah initio MD has mostly relied on two 
general approaches: Born-Oppenheimer MD (BOMD) 
and Car-Parrinello MD (CPMD), each with its advan- 
tages and shortcomings. In BOMD, the electronic wave 
function energy is iteratively minimized for each MD time 
step, making this approach computationally rather ex- 
pensive. The CPMD&i£ approach bypasses the expen- 
sive iterative minimization by introducing the fictitious 
mass for the electronic degrees of freedom and a suit- 
ably designed "on-the-fly" scheme for their propagation. 
The fictitious mass has to be small enough to ensure that 
the electrons follow the nuclei adiabatically, very close to 
their instantaneous electronic ground stated. As a conse- 
quence, the maximum permissible integration time step 
in CPMD is significantly smaller than that of BOMD, 
thus limiting the simulation timescales. 

The recently developed second-generation CPMD 
method combines the best of both approaches by retain- 
ing the large integration time steps of BOMD and, at the 
same time, preserving the efficiency of CPMDi 4 -. The 
propagation scheme in second-generation CPMD relies 
on the appropriate usage of the information about the 
electrons from the previous MD steps and, in contrast 
to CPMD, does not require a fictitious mass parameter 
to maintain the accuracy of the Born-Oppenheimer dy- 
namics. The superior efficiency of this new approach, 
which, depending on the system, varies between one to 
two orders of magnitude, has been demonstrated for a 
large variety of different systems^—. 

Propagation of the electronic degrees of free- 
dom. Within mean-field electronic structure methods, 
such as Hartrec-Fock and Kohn-Sham density functional 
theory (DFT), the electronic wave function is described 



by a set of occupied molecular orbitals (MOs) or 
by an associated idempotent one-electron density oper- 
ator p = J^i In the second-generation CPMD 
method, the propagation of the electronic degrees of free- 
dom is achieved by adapting the predictor-corrector inte- 
grator of Kolaf a 21 ' 22 to the electronic structure problem. 
Firstly, the predicted MOs at time t n are constructed in 
terms of the electronic degrees of freedom from the K 
previous MD steps: 
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Secondly, the orbitals are corrected by performing a sin- 
gle step \5ip p (t n )) along the preconditioned electronic 
gradient direction computed with the orbital transforma- 
tion method 2 ^. This leads to the final updated orbitals: 
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The second-generation CPMD method avoids the self- 
consistency cycle entirely and obviates the need for 
the computational expensive diagonalization, replacing 
it with a fast orbital transformation step. The electron 
propagation scheme is rather accurate and time reversible 
up to 0(At 2K ~ 2 ), where At is the integration time step. 
It keeps the electronic degrees of freedom very close to 
the Born-Oppenheimer values and allows for At to be as 
large as in BOMD. 

Energy functional and nuclear forces. The to- 
tal energy E PC [p p ] is evaluated from both predicted and 
corrected orbitals and represents an approximation to the 
Harris-Foulkes functiona l 24 ' 25 : 

,P p (r)p p (r') 
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where p p {r) is the predicted electron density associated 
with \ip p (t n )), H[p p ] is the effective Hamiltonian oper- 
ator, vxc[p p ] the exchange-correlation (XC) potential, 
whereas Exc[p p ] an d En are the XC energy and the 
nuclear Coulomb interaction, respectively. 

The nuclear forces are computed by evaluating the an- 
alytic energy gradient Ff c = Vr i .E pc '[/> p ]. However, 
since Ap = p — p p ^ 0, the usual Hellmann-Feynman 2 ^ 
and Pulay forces 2 ^ have to be augmented by the following 
extra term: 
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where Vh[p p ] is the Hartrcc potential, while p is the cor- 
rected density that corresponds to \tpi(t n )). 
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Modified Langevin equation. Despite the close 
proximity of the electronic degrees to the Born- 
Oppcnhcimcr surface, the employed propagation scheme 
introduces a small dissipative energy drift during long 
MD runs. However, it is possible to correct for that by 
devising a modified Langevin-type equation that in its 
general form reads as: 

MtRj = Ff° - jMiRi + Si, (5) 

where Mi are the nuclear masses, Ri the nuclear coor- 
dinates (the dot denotes time derivative), Ff° the ex- 
act Born-Oppenheimer forces, while 7 a damping coef- 
ficient, and H/ is an additive white noise, which must 
obey the fluctuation-dissipation theorem (Ej(O)Ej(t)) = 
2 n fkbTMiS(i) in order to sample the canonical distribu- 
tion. 

Presuming that the energy is exponentially decay- 
ing, which had been shown to be an excellent assump- 
tio n 12 ' 14 ' 28 , it is possible to rigorously correct for the dis- 
sipation, by modelling the nuclear forces arising from our 
dynamics as: 

Ff c = Ff° - lD MiRi, (6) 

where jd is an intrinsic, yet unknown damping coefficient 
to mimic the dissipation. 

By substituting Eq. (|6|) into Eq. ([5]) the following mod- 
ified Langevin-like equation is recovered: 

M1R1 = Ff c + H/ (7) 

In other words, if we would know the unknown value of 
7b, in spite of the dissipation, it is nevertheless possible 
to guarantee an exact canonical sampling of the Boltz- 
mann distribution, by simply augmenting Eq. [5] with E / 
according to the fluctuation-dissipation theorem. For- 
tunately, the intrinsic value of 7^ does not need to be 
known a priori, but inspired by ideas of Krajewski and 
Parrinello^, can be determined in a preliminary run us- 
ing a Berendsen-like algorithm^ in such a way that even- 
tually the equipartition theorem (^MiRj\ = \kpT 
holds. Although this can be somewhat laborious, but 
once 7d is determined, very long and accurate simula- 
tions can be routinely performed at a greatly reduced 
computational cost. 

B. Decomposition analysis based on absolutely 
localized molecular orbitals 

Energy decomposition analysis. One of the most 
powerful techniques that modern first-principle electronic 
structure methods provide to study and analyse the in- 
termolecular bonding is the decomposition of the to- 
tal molecular binding energy into physically meaningful 
components. Such methods have shown that the HB 
is a result of the delicate interplay between weak dis- 
persive forces, electrostatic effects (e.g. charge-charge, 
charge-dipole and charge-induced dipole interactions) 



and donor-acceptor type orbital (i.e. covalcnt) interac- 
tions, such as forward- and back-donation of electron 
density between water molecules. The extent of the dif- 
ferent modes of interactions determines the strength of 
the HBs in water clusters and condensed phases and, con- 
sequently, all their physical properties. 

The need for physically reasonable and quantitative 
useful values of the energy components has resulted in 
numerous decomposition schemes which have been pro- 
posed since the early years of theoretical chemistry2Ir— . 
In this review, we discuss a novel energy decomposition 
analysis scheme based on absolutely localized molecular 
orbitals (ALMO EDA) 4 ^. Unlike conventional MOs, 
which are generally delocalized over all molecules in the 
system, ALMOs are expanded in terms of the atomic 
orbitals of only a given molecule^— . Although AL- 
MOs were originally used to speed up the calculation 
of SCF energies for large ensembles of molecules^, they 
arc now widely used in energy decomposition analysis 
(EDA)iiri£. It should be mentioned that since the in- 
troduction of ALMO EDA4£, ALMO-based decomposi- 
tion methods have been extended to many-determinant 
wave functions^. Nevertheless, the application of the 
most recent approach to water is still very limited and 
here we focus on the decomposition procedure based on 
the mean-field methods such as Hartree-Fock and density 
functional theory. 

ALMO EDA separates the total interaction energy 
of molecules {AEtot) into a frozen density component 
(FRZ), as well as polarization (POL) and electron dclo- 
calization (DEL) terms 

AEtot = AEfrz + AEpol + AEdel- (8) 

The frozen density term is defined as the energy 
change that corresponds to bringing infinitely separated 
molecules into the complex geometry without any relax- 
ation of the MOs on the monomers, apart from energetic 
modifications associated with satisfying the Pauli exclu- 
sion principle: 

AE FRZ = E(RpKz) - ]T E(R X ), (9) 

X 

where E(R X ) is the energy of the isolated molecule x with 
its nuclei fixed at the geometry that it has in the system. 
Rfrz is the density matrix of the system constructed 
from the unrclaxcd MOs of the isolated molecules. 

The polarization energy is defined as the energy 
lowering due to the intramolecular relaxation of each 
molecule's ALMOs in the field of all other molecules in 
the system. The intramolecular relaxation is constrained 
to include only variations that keep MOs localized on 
their molecules, i.e. 

AEpol = E(Rpol) - E(Rfrz), (10) 

where Rpoi is the density matrix constructed from the 
fully optimized (polarized) ALMOs. That is, the Rpol 
state is the lowest energy state that can be constructed 
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from completely localized MOs. ALMOs are not orthog- 
onal from one molecule to the next and, therefore, the 
minimization of the electronic energy as a function of the 
ALMO coefficients differs from conventional SCF meth- 
ods. Mathematical and algorithmic details of the self- 
consistent field procedure for finding the polarized state 
Rpol have been described by many authors 4 ^— . 

The remaining portion of the total interaction en- 
ergy, the electron derealization (DEL) or charge-transfer 
(CT) energy term, is calculated as the energy difference 
between the state formed from the polarized ALMOs 
{Rpol) and the state constructed from the fully opti- 
mized delocalized MOs (Rscf)- Thus, 

AEdel = E{R SC f) - E{R POL ), (11) 

where AEdel includes the energy lowering due to elec- 
tron transfer from the occupied ALMOs on one molecule 
to the virtual orbitals of another molecule, as well as the 
further energy change caused by induction (or repolariza- 
tion) that accompanies such an occupied-virtual orbital 
mixing. Therefore, AEdel is further decomposed into 
the occupied- virtual forward- and back-donation compo- 
nents for each pair of molecules, plus the higher order 
(HO) induction energy: 

AEdel = i AE ^v + AE v^} + ae ho, (12) 

x,y<x 

AEho is generally small and cannot be naturally divided 
into two-body terms. 

Charge transfer analysis. In addition to the en- 
ergy decomposition, ALMOs have been used to mea- 
sure the amount of electron density transferred between 
molecules^. The intermolecular CT in such a charge 
transfer analysis (ALMO CTA), AQdel, is defined as 
the change in the electron density from the polarized 
Rpol state to the fully converged state Rscf- This defi- 
nition is well-justified because, according to the Mullikcn 
analysis, the ALMO constraint explicitly excludes CT 
between molecules, making Rpol the "zero-CT" state 
with the lowest energy. 

In agreement with ALMO EDA, AQdel includes the 
charge transfer due to occupied-virtual mixing (AQx^y) 
and the accompanying higher order relaxation terms 
(AQho)- 

AQdel = J2 i A Q^y + A Qy^} + a Qho- (13) 

x : yKx 

Thus, the intermolecular charge transfer terms in ALMO 
CTA have corresponding well-defined energies of stabi- 
lization calculated by ALMO EDA. 

Significant complementary occupied-virtual 
pairs. In addition to quantifying the amount and 
energetics of intermolecular charge transfer, it is often 
useful to have a simple description of orbital interactions 
between molecules. The polarized canonical ALMOs 
used as a reference basis set in the decomposition anal- 
ysis do not directly show which occupied- virtual orbital 



pairs are the most important in forming intermolecular 
bonds. That is, in general there are no occupied-virtual 
pairs in this basis set that can be neglected. However, 
orbital rotations within the occupied subset and within 
the virtual subset of a molecule leave AE x ^ y and 
AQ x ^y unchanged. This freedom can be used to find 
new sets of orbitals for x and y, in which charge transfer 
from x to y is described as each occupied orbital on x 
donating electrons to only one (complementary) virtual 
orbital on y. Such orbitals are called complementary 
occupied-virtual pairs (COVPs)^. 

Construction of the COVPs greatly simplifies the pic- 
ture of intermolecular orbital interactions since they form 
a "chemist's basis set" for a conceptual description of 
bonding in terms of just a few localized orbitals. As we 
demonstrate below, COVPs provide an alternative and 
somewhat unconventional view of hydrogen bonding in 
the water dimer. 

Computational details and implementation. 

ALMO EDA and CTA are implemented in both, the Q- 
Chern^ and CP2 K 53 i 54 software packages. The varia- 
tional optimization of the occupied ALMOs is performed 
using the locally projected SCF metho d 48 i 50 . The defi- 
nition of the ALMOs and the polarization energy relies 
on an underlying basis set that is partitioned amongst 
the fragments. Gaussian AO basis sets in both packages 
are ideal in this regard, and give well-defined polarization 
energies as long as there are no linear dependences. In 
the linearly dependent limit, where the basis functions 
on one fragment can exactly mimic functions on another 
fragment, this ceases to be the case^i. This is not an 
issue for the AO basis sets used routinely in quantum 
chemistry and AIMD simulations^. 

In addition to the occupied ALMOs, the locally pro- 
jected SCF method yields a set of non-redundant linearly 
independent virtual ALMOs. After the SCF procedure is 
converged, the occupied subspacc is projected out from 
the virtual ALMOs to ensure strong orthogonality of the 
subspaces. The energy lowering due to electron transfer 
from the occupied ALMOs on molecule x to the virtual 
orbitals of molecule y in Eq. [T^] is a quasi-perturbative 
energy correctio n 50 i 56 i 57 that can be expressed as: 

Ox V V 

AE x ^ y = J2J2 FXi yaX ya xi , (14) 

i a 

where F xl ya is the contravariant-covariant representation 
of the Fock operator build from the converged ALMOs 
and X ya X i is the amplitude corresponding to the elec- 
tron transfer (excitation) from the converged absolutely 
localized occupied orbital i on fragment x to the virtual 
orbital a on fragment y. The variational nature of the po- 
larized ALMOs guarantees that the energy term within 
a molecule is zero, AE X ^ X = 0. 

The corresponding amount of charged transferred from 
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x to y (Eq. [T3|) is expressed as^ 6 -: 

Ox V V 

AQ^ y = E A'",,,.V"-,, (15) 

z a 

The amplitudes X^ ^ are obtained by solving a 
quadratic equation^ 7 - using the preconditioned conjugate 
gradient method: 

V o 

J?Va . i \ " pya vwb . _ \ " yya pzj 
O V 

EE Xy a z] F^ wb X wb xl = (16) 

zj wb 

The total energy lowering due to the occupied-virtual 
mixing J2 X AE x ^ y is equivalent to the result obtained 
by a single Fock matrix diagonalizatio n 50 i 56 . The AEho 
term is introduced to recover the small difference between 
this energy and the fully converged SCF result. 

It is important to note that special care must be 
taken to remove the basis set superposition error (BSSE) 
from the interaction energies and their components. The 
BSSE is not introduced when calculating frozen den- 
sity and polarization energy contributions because con- 
strained ALMO optimization prevents electrons on one 
molecule from borrowing the AOs of other molecules to 
compensate for the incompleteness of their own AOs. 
However, the BSSE enters the charge transfer terms 
since both the BSSE and charge transfer results from 
the same physical phenomenon of derealization of frag- 
ment MOs. Therefore, these terms are inseparable from 
each other when finite Gaussian basis sets are used to 
describe fragments at finite spatial separation. It has 
been demonstrated that the BSSE decreases faster than 
charge transfer effects with increasing quality of the ba- 
sis S e t 50 i 58 i 59 . Therefore, the use of medium and large lo- 
calized Gaussian basis sets (without linear dependencies) 
make the BSSE component of the interaction energy neg- 
ligibly small. The BSSE associated with each forward- 
and back-donation term AE x ^, y can be corrected indi- 
vidually, as shown in Refs. HH andHfl 

Features of the ALMO decomposition methods. 
ALMO EDA is conceptually similar to long-established 
decomposition methods, such as the Morokuma analy- 
sis^ 3 -, but includes several important novel features de- 
scribed below. 

• Unlike earlier decomposition methods 3 ^—, ALMO 
EDA and CTA treat the polarization term in a 
variationally optimal way. Therefore, CT effects 
(i.e. effects due to intermolecular electron dereal- 
ization) cannot be over- or underestimated. 

• The CT term can be decomposed into forward- 
donation and back-bonding contributions for each 
pair of molecules in the systems. 



• The ALMO charge transfer scale is such that all 
terms have well defined energetic effects. In con- 
trast, population analysis methods include not only 
the true CT, but also large contaminating charge 
overlap effects^ 6 -. 

• COVPs constructed from canonical ALMOs pro- 
vide a compact and chemically intuitive description 
of electron transfer between the molecules. 

• The ALMO method in the CP2K package is cur- 
rently the only decomposition scheme for con- 
densed matter systems. CP2K relies on the mixed 
Gaussian and plane wave representation of elec- 
trons^, which makes it uniquely suited for per- 
forming ALMO calculations for periodic systems. 
In CP2K, the localized atom-centered Gaussian ba- 
sis sets are used for the construction of ALMOs, 
whereas plane waves ensure the computational effi- 
ciency in large-scale calculations of the Hartree and 
XC potentials for periodic systems. 



III. PHYSICAL NATURE OF HYDROGEN 
BONDING IN THE WATER DIMER 

As already alluded to, hydrogen bonding is central 
to all aqueous systems ranging from water nanoclus- 
ters and microsolvated ions to bulk water and solvated 
biomolccules^— . Despite numerous experimental and 
theoretical studies 6 ^—, the physical nature of hydrogen 
bonding is still under debate. One issue is the degree of 
covalency in the hydrogen bonding, which is determined 
by the extent of intermolecular electron derealization 
or CTlSrll Natural bond orbital (NBO)S and natural 
EDA2 8 . suggests that CT is predominan t 39 ' 40 ' 76 because 
if CT is neglected, NBO analysis shows no binding at 
the water dimer equilibrium geometry. However, other 
earlier decomposition methods 3 ^— estimated that CT 
contributes only around 20% of the overall binding en- 



The extent of CT has practical significance for aque- 
ous molecular dynamics simulations, where models based 
on purely electrostatic potentials (e.g. Coulomb per- 
haps with polarizability plus Lcnnard- Jones) seem to be 
very successful in reproducing many of the structural 
and thermodynamic properties of wate r 78 ' 79 . However, 
the failure of classical molecular dynamics simulations 
to reproduce the controversial "chains and rings" struc- 
ture of liquid water, inferred from recent X-ray absorp- 
tion and X-ray Raman scattering experiments^, has gen- 
erated questions about the reliability of existing water 
potcntial o 2 ' 76 ' 80 . This fact, combined with the predomi- 
nantly CT character of the hydrogen bonding suggested 
by NBO, has led to proposals to incorporate CT effects 
into empirical water potentials 7 - 6 -. 

In this section, we review the role of intermolecular CT 
effects in the simplest water cluster - the water dimer - 
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uncovered with ALMO EDA and CTA^i. Accurate sepa- 
ration of polarization effects from CT is essential to deter- 
mine the amount of covalency in the hydrogen bonding. 
A variationally optimal treatment of polarization makes 
ALMO EDA and CTA ideal for this purpose. 

Energetic components of the hydrogen bond 
stabilization. The relative position of the molecules 
in the water dimer with C s symmetry is described by 
three parameters shown in Figure [T]A The structure 
of the dimer was optimized at the MP2/aug-cc-pVQZ 
level and is characterized by a = 172°, = 127°, and 
Roh = 1.94A As mentioned above, the ALMO decom- 
position analysis is presently limited to single determi- 
nant wave functions. Therefore, ALMO EDA was ap- 
plied to wave functions calculated at the Hartree-Fock 
and DFT level using a series of local density approxima- 
tion (LDA), generalized gradient approximation (GGA), 
hybrid and meta- hybrid XC functionals (Table HJ. 

The decomposition of the Hartree-Fock energy pro- 
duces results similar to the earlier decomposition meth- 
ods, but gives a somewhat larger charge transfer con- 
tribution (Table HJ). According to ALMO EDA, charge 
transfer accounts for 27% of the total Hartree-Fock bind- 
ing energy. When Kohn-Sham DFT is used instead of the 
Hartree-Fock method (Table IJ), all energy terms change 
because of modification of the exchange and addition of 
the correlation terms into the mean-field Hamiltonian. 
The derealization effect becomes more pronounced for 
the density functional methods and in some cases the 
charge transfer term is more than 45% of the overall bind- 
ing energy. This observation is consistent with the ten- 
dency of modern density functionals to underestimate the 
HOMO-LUMO gap 82 i 83 , which in the water dimer case, 
manifests itself in a larger charge transfer energy. 

It is clear from Table fl] that all three energy compo- 
nents (frozen density, polarization, and charge-transfer) 
are important for the energetic stabilization of the dimer 
at its equilibrium geometry. We, therefore, conclude that 
the NBO approach significantly overestimates CT due to 
the non- variational treatment of the reference "zero-CT" 
electronic state. Our results show that CT contributes 
around one third of the overall binding energy in the 
complete basis set limit, of which approximately 95% is 
from the proton acceptor to the proton donor. The same 



FIG. 2. Dependence of the energy components on the dis- 
tance between the water molecules in the water dimer at the 
HF/aug-cc-pVQZ level of theory. 



effect is observed on the charge scale, indicating a direct 
correspondence between electron redistribution and the 
energy of CT interactions. 

The relative contribution of the energy terms varies 
strongly with the position of the molecules in the water 
dimer. Figure [2] shows the dependence of the Hartree- 
Fock BSSE corrected energy and its ALMO decomposi- 
tion on the distance between the water molecules (Ron 
is varied and all other internal coordinates remain fixed 
at their MP2/aug-cc-pVQZ values). The frozen density 
component increases significantly and becomes repulsive 
as the molecules get closer. At the same time, stabi- 
lization due to polarization and charge transfer increases 
upon the closer contact, but not strongly enough to com- 
pensate for the electron density repulsion. The stabiliz- 
ing contribution of charge transfer and polarization de- 
crease rapidly with the increase of the intermolecular dis- 
tance and, at Roh > 3A, the interaction energy can 
be accurately approximated by the frozen density term 
alone. 

Table|T]shows that the relative contribution of the three 
energy components somewhat depends on the XC func- 
tional chosen, just as the total binding energies do. Such 
a wide spread of the DFT results indicates that a bet- 
ter description of the intermolecular correlation energy 
is necessary for the water dimer. At the same time, 
the main qualitative description of hydrogen bonding re- 
mains the same for all commonly used XC functionals. 
Therefore, we further discuss the results obtained with 
the B3LYP XC functional. This functional most closely 
reproduces more accurate MP2 water dimer binding en- 
ergies of the various XC functionals that we tried (Ta- 
ble IJ) . A detailed discussion of the performance of popu- 
lar density functionals for describing HB between water 
molecules can be found in Refs. 1921494 

The results of the B3LYP ALMO decomposition analy- 
sis are presented in Table [TH which shows that the energy 
and charge components rapidly converge as the basis set 
becomes locally complete, indicating the high stability of 
the ALMO decomposition. All CT terms presented here 
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TABLE I. The decomposition analysis results are reported for an LDA (SVWN^M^), three GGA (BPSr^^ 7 ., PW91^), a hybrid 
and an meta-hybrid XC functional (BMK—), as well as for the Hartree-Fock wave function. All terms are BSSE 
corrected and an aug-cc-pVQZ basis set is employed, throughout. 



Scale AQ, me AE, kj/mol 



Functional 


HF 


BP86 


BMK 


B3LYP 


PW91 


SVWN 


HF 


BP86 


BMK 


B3LYP 


PW91 


SVWN 


FRZ 








0.0 






-5.0 


-1.8 


-5.4 


-5.2 


-8.5 


-16.6 


POL 








0.0 






-6.1 


-7.0 


-6.6 


-6.5 


-5.2 


-7.9 


A^D a 


0.1 


0.2 


0.1 


0.1 


0.2 


0.2 


-0.2 


-0.3 


-0.2 


-0.3 


-0.4 


-0.3 


d^aI 


0.8 


4.0 


1.4 


2.8 


4.7 


4.3 


-3.2 


-8.1 


-5.1 


-6.6 


-8.1 


-8.1 


Rem.~CT^ 


0.2 


-1.1 


0.3 


-0.4 


-1.8 


-1.4 


-0.6 


0.1 


-0.7 


-0.3 


0.1 


0.1 


TOT_ 


1.1 


3.1 


1.8 


2.5 


3.1 


3.1 


-15.1 


-17.1 


-17.9 


-18.9 


-22.1 


-32.8 


bsse" 


0.0 


0.0 


0.1 


0.0 


0.0 


0.0 


0.1 


0.2 


0.6 


0.1 


0.2 


0.2 


COVPi^. 


06 


07 


83 


96 


95 


96 


94 


95 


85 


93 


90 


93 



a D - electron-donor (proton-acceptor), A - electron- acceptor (proton-donor). 

b Remaining CT includes intramolecular terms as well as the higher order relaxation term. 

c BSSE corrected MP2/aug-cc-pVQZ total interaction energy is -20.6 kj/mol. 

d Contribution of COVPi is given as percent of D — > A. 



are BSSE corrected. The BSSE is presented in Table ED 
to show the degree of basis set completeness. The very 
small values of BSSE suggests that the aug-cc-pVQZ and 
aug-cc-pV5Z basis sets are effectively complete for the 
present purpose. 

The role of electron transfer in hydrogen bond- 
ing. The total electron density transfer calculated with 
ALMO CTA is just a few milli-electrons (2.3me at the 
B3LYP/aug-cc-pV5Z level). This result is an order 
of magnitude smaller than the CT calculated with the 
Mullikcn, Ldwdin and natural population analysis (PA) 
methods (Table . This discrepancy arises largely 

from the different meaning assigned to CT in ALMO 
CTA and PA techniques. ALMO CTA measures CT as 
the degree of electron relaxation from the optimal polar- 
ized (pre-CT) state to the delocalized state. By contrast, 
PA methods include not only the true CT, but also the 
separate and, in this case, larger effect of partitioning the 
charge distribution of the polarized pre-CT state (for a 
detailed comparison see Ref. [iih . Thus, the key advan- 
tage of the ALMO CTA approach is that it shows the 
electron transfer associated with an energy lowering due 
to dative interactions: just a few milli-electrons. 

It may seem remarkable that so little CT can stabilize 
the HB by 6.5 kj/mol (equivalent to 32 eV per electron), 
but this estimate is consistent with simple estimates from 
perturbation theory. The CT energy is a second order 
correction to the energy of the polarized system, and is 

F 2 

proportional to -, — ^-^r , where F a d is the CT energy cou- 
pling between donating orbital d and accepting orbital a, 
while €i is the energy of orbital i. The AQ term, however, 

F 2 

is proportional to , _ ad -, 2 . Therefore, the CT energy per 
electron is related to the energy gap between donating 
and accepting orbitals e a — e<i- The B3LYP/aug-cc-pV5Z 
energy gap between the most important donating and ac- 
cepting orbitals in the water dimer lies between 10 and 
40 eV (virtual orbitals in such a big basis set practically 
form a continuum of states). Thus, a value of 32 eV for 



= 60° = 127° 

'HQ iH-Qfr 

= 179° = 240° 

FIG. 3. Dependence of the shape of the most significant 
COVP on the relative orientation of the water molecules in the 
dimer. All figures show orbitals calculated at the B3LYP/aug- 
cc-pVQZ level of theory. Isosurface value of 0.05 a.u. Occu- 
pied orbitals are represented with saturated colours. Faint 
colours represent complementary virtual orbitals. 



the effective d — a gap for CT between water molecules 
in the dimer at the equilibrium geometry is reasonable. 

Orbital representation of donor-acceptor inter- 
actions. COVPs let us visualize CT effects and provide 
additional insight into the nature of hydrogen bonding. 
In the water dimer, only one COVP is significant, and 
recovers 96-97% of the overall transfer an from proton ac- 
ceptor to a proton donor on the energy and charge scales 
(Table [II]). The remaining CT in this direction can be 
attributed to the four remaining COVPs, none of which 
exceeds 3% of the overall transfer. The shapes of the 
occupied and the virtual orbitals of the main COVP are 
shown in Figure [3] (O = 127°). The virtual orbital resem- 
bles the O-H anti-bonding orbital, <J* OHl of the electron 
accepting molecule. This is consistent with oxygen atom 
lone pairs donating electron density to the anti-bonding 
orbitals of the other molecule. However, the shape of the 
donating orbital is somewhat unexpected because it does 
not resemble an sp 3 -hybridiscd lone pair. 
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TABLE II. ALMO CTA and EDA results for the water dimer (B3LYP/aug-cc-pVYZ). All terms are BSSE corrected. 



Scale 






AQ, me 








AE, kJ/mol 




Y 


D 


T 




Q 


5 


D 


T 


Q 


5 


FRZ 






0.0 






-5.5 


-5.4 


-5.2 


-5.1 


POL 






0.0 






-4.5 


-6.1 


-6.5 


-7.1 


A^D^ 


0.1 


0.2 




0.1 


0.1 


-0.2 


-0.4 


-0.3 


-0.2 


D^A^ 


3.7 


2.4 




2.8 


2.7 


-7.9 


-6.6 


-6.6 


-6.3 


Rem."c r f* 


0.3 


0.2 




-0.4 


-0.5 


-0.4 


-0.3 


-0.3 


-0.3 


Tor 


4.0 


2.7 




2.5 


2.3 


-18.5 


-18.8 


-18.9 


-18.9 


bsse" 


0.3 


0.1 




0.0 


0.0 


1.0 


0.2 


0.1 


0.0 


COVP^ 


95 


97 




96 


97 


90 


97 


93 


96 



a D - electron-donor (proton-acceptor), A - electron- acceptor (proton-donor). 
b Remaining CT includes intramolecular terms as well as the higher order relaxation term. 
c MP2 interaction energies are -18.3, -19.8, and -20.6 kJ/mol for Y=D, T, Q, respectively. 
d Contribution of COVPi is given as percent of D — > A. 



TABLE III. Charge (me) of the electron-acceptor molecule 
in the dimer (B3LYP/aug-cc-pVXZ). All charges are BSSE 
corrected. 



X 


D 


T 


Q 


5 


Mulliken PA 


-27.5 


-18.8 


-22.2 


-17.0 


Lowdin PA 


-24.0 


-24.0 


-21.2 


-17.8 


Natural PA 


-18.3 


-16.5 


-17.1 





Commonly, sp 3 -hybrids play an important role in pre- 
dicting geometries of gas-phase water molecules and di- 
rection of HBs in ice phases. Furthermore, sp 3 -hybridised 
lone pairs on the O atom have become a commonly ac- 
cepted way to visualize the electronic structure of water 
molecules and are the basis for the five-point molecular 
models (e.g. TIP5P) widely used in classical molecular 
dynamics simulations^. However, the occupied orbitals 
are not unique, and can, therefore, be fixed by criteria 
that include well-defined ionization energy (giving canon- 
ical orbitals, Figure [4]), maximal localization (giving sp 3 
lone pairs), or the most compact representation of CT ef- 
fects (giving the COVP shown in Figure The COVP 
is best suited for studying donor-acceptor interactions, 
and the form of the optimal donor and acceptor orbitals 
can be understood as a compromise between high energy 
and good interaction with the acceptor. In this regard, 
the optimal acceptor orbital of Figure [3] bears almost no 
resemblance to the low-lying canonical virtual orbitals of 
Figure SI consistent with the effective d — a gap being far 
larger than the HOMO-LUMO gap. The occupied (do- 
nating) orbital is mostly a linear combination of 3ai and 
l&i canonical orbitals (Figure that are the two highest 
lying orbitals of the electron donating water molecule. 
A small CT contribution from the 2a\ canonical orbital 
of the donating molecule is reasonable given that its low 
energy makes it a poor donor. This simple argument ex- 
plains the form of the donating orbital in the water dimer 
complex (Figure El 8 = 127°). 

Further support for this interpretation comes from Fig- 
ure [SJ which shows the dependence of the B3LYP/aug-cc- 



1a 1 (-520.7) 23,1-27.7) 1b 2 (-14.7) 3a, (-10.9) 1b, (-8.8) 



A A 

4a, (-0.7) 2b 2 (0.5) 5a, (1.7) 



2b, (2.2) 



FIG. 4. Symmetry of the five occupied and the lowest four 
virtual canonical MOs of water molecule. See the caption of 
Figure [3] for full description. 



pVQZ BSSE corrected energy and its ALMO decomposi- 
tion on the orientation of the water molecules (0 is varied 
and all other internal coordinates remain fixed at their 
MP2/aug-cc-pVQZ values). The CT energy does not 
maximize around tetrahedral coordination (0 = 127°). 
The energy lowering due to the most significant COVP 
changes remarkably little from -4.9 kJ/mol for O = 180° 
to -7.7 kJ/mol for 8 = 50°. From Figure|3j the donor or- 
bital does not rotate with the water molecule but stays di- 
rected towards the electron accepting molecule, unlike an 
sp 3 lone pair. The principal donor orbital thus changes 
with rotation to optimize the coupling with the comple- 
mentary ctq H virtual acceptor orbital, thereby explaining 
the weak dependence of the CT energy on O. 

It is interesting to consider HBs that involve bifurcated 
interactions. O = 179° corresponds to an OH bond inter- 
acting with two sp 3 lone pairs, which represent a bifur- 
cated proton donor in the traditional picture (see the car- 
toon in Figure QJ3). However, the CT contribution to H- 
bonding still involves only one donating orbital, as shown 
in Figure [3] and the cartoon in Figure [TC. The reduction 
in CT energy reflects at O = 179° a greater contribution 
of the lower energy 3ai orbital and a decreased contri- 
bution of the highest occupied lbi orbital of the donor 
molecule. In fact, the CT energy dependence on O can 



E -5.0 
3 



-1>FRZ + P0L 
o CTduetoCOVPl 




110 130 150 



170 190 





210 230 250 270 290 



FIG. 5. Dependence of the B3LYP/aug-cc-pVQZ energy com- 
ponents on the relative orientation of water molecules in the 
dimer. Dashed line represents CT calculated according to 
Eq. 1171 Dotted line represents FRZ+POL interactions calcu- 
lated according to Eq. [18] 



be well represented as a linear combination of CT from 
these two orbitals by a simple equation (the dashed line 
in Figure [5]): 

AE D ^ A (e) = AE lbl sin 2 (9) + AE 3ai cos 2 (9) (17) 

The shape of the FRZ curve can be explained in purely 
electrostatic terms. The dotted line in Figure [5] repre- 
sents the interaction energy of point charges placed at 
the position of the nuclei in the dimer (— l.Oe and +0.55 
charges replace C and H atoms correspondingly) and is 
given by equation: 

AE FRZ+POL (Q) = ^^2^+ 29 ' 2 kJ m ° rl ( 18 ) 

The constant in the equation is included to capture 
effects that are essentially independent of 9 such as po- 
larization and exchange. Therefore, the position of the 
minimum on the total energy curve at 9 = 127° is deter- 
mined by combination of both electrostatic and charge 
transfer interactions. 

In the case of a bifurcated proton acceptor (Figure[Tp), 
the description of HB changes qualitatively. The CT 
term becomes very small due to poor interactions and 
has significant contributions from two CCVPs - the sym- 
metric (3ai canonical orbital) and anti-symmetric (l&i 
canonical orbital) (Figures @] and . This reflects the 
availability of two acceptor anti-bonding o* OH (symmet- 
ric 4ai and anti-symmetric 262) orbitals in the vicinity 
of the electron donating molecule. It is certain that the 
donating orbitals will change their shape and orientation 
in larger water clusters and bulk liquid water according 
to the local environment. 

Summary. ALMC energy decomposition and charge 
transfer analysis is a promising method to study hydro- 
gen bonding between water molecules using density func- 



COVPi 
64% of AE CT{D ^ A) 
59% of AQct (D -» a} 



COVP 2 
29% of AE CT(D ^ A) 
37% of AQa{D-»A) 



FIG. 6. Bifurcated hydrogen bonding in the water dimer. 
EDA terms, kj/mol: AE FRZ = -8.2, AE PO l = -1.7, 
AEd^a = -0.9, AEa^d = -0.2, AEtot = -11.0. CTA 
terms, me: AQ D ^ A = 0.4, AQa^d = 0.1, AQtot = 0.5. 
See the caption of Figure [3] for full description. 



tional theory. It shows that although electron derealiza- 
tion effects play an important role in hydrogen bonding 
they are not solely responsible for the energetic stabiliza- 
tion of the water dimer. The contributions of frozen den- 
sity interactions and polarization are not less significant 
than that of CT, unlike some earlier work^^^. ALMC 
CTA demonstrates that the amount of intcrmolccular 
CT is in the order of few milli-electrons, which is much 
smaller than has been inferred from population analyses. 
It also shows that CT is fairly insensitive to intermolecu- 
lar rotation of the water molecules. This helps to account 
for the success of empirical potentials that do not include 
charge transfer explicitly. Furthermore, COVPs provide 
a new view of the electron donating orbital in the wa- 
ter dimer. Unlike rigid sp 3 lone pairs, the COVP donor 
changes its orientation according to the relative positions 
of the two molecules. A single p-likc lone pair is usually 
the dominant donor, although at the geometry of a bifur- 
cated HB, the CT contribution becomes small and two 
donor orbitals contribute. 

Further applications of ALMO EDA to water 
clusters. ALMC EDA has been applied to investigate 
the importance of charge transfer effects for the vibra- 
tional spectrum of the water dimer—. Comparing the vi- 
brational spectra calculated with and without the charge- 
transfer shows that electron derealization has a very 
large effect on the vibrational frequency and intensity 
associated with the stretch of the donated OH bond, 
while its effect on the other vibrational modes is small2&. 
Further applications of ALMO EDA to water clusters 
include a recent study of HB in trimer, tetramer, pen- 
tamer, 13-er, and 17-er— . The decomposition of the 
two- three- and higher-body interaction energies into the 
frozen-density, polarization and charge-transfer compo- 
nents has revealed several interesting trends that provide 
additional physical justification for the standard practice 
of not explicitly including charge transfer into water force 
fields. ALMO EDA has also been used, in combination 
with other decomposition schemes, to gain insight into 
the performance of several popular density functionals 
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for describing interactions between water molecules^ 3 -. 



IV. LIQUID WATER 

A. Structural and dynamical properties of liquid 
water 

Despite the ongoing development of new simulation 
techniques, an accurate modelling of liquid water still 
represents a major challenge. The accuracy of AIMD 
simulations has often to be sacrificed to reduce their com- 
putational burden. Therefore, various physical effects 
such as van der Waals interactions between molecules 
and the quantum behaviour of nuclei have to been repro- 
duced only approximately or even completely neglected. 
In addition, the high computational cost of AIMD im- 
poses severe restrictions on the size of a model system 
and time length of the simulations thus introducing ad- 
ditional finite-size errors and statistical uncertainties into 
the properties calculated^r-iSi. 

In this section, we assess the magnitude of errors 
coming from physical approximations, finite-size effects 
and insufficient sampling by performing large-scale sim- 
ulations, which exploit the computationally efficiency 
second-generation CPMD metho d 14 ' 108 . 

Computational details. The largest simulated sys- 
tem consisted of 128 light water molecules in a periodic 
cubic box of length L = 15.6627 A, which corresponds to 
a density that differs by only 0.3 % from the experimental 
value. All simulations were performed at 300 K within 
the canonical NVT ensemble; the Langevin equation of 
motion was integrated using the algorithm of Ricci and 
Ciccottii ^. The discretized integration time step At was 
set to 0.5 fs, while j D = 8.65 x 10" 5 fs" 1 . The simul- 
taneous propagation of the electronic degrees of freedom 
proceeded with K = 7, which yields a time reversibility 
of 0(At 12 ). At each MD step the corrector was applied 
only once, which implies just one preconditioned gradi- 
ent calculation. The deviation from the BO surface, as 
measured by the preconditioned mean gradient deviation 
was 10" J a.M., which is only slightly larger than typical 
values used in fully converged BOMD simulations. The 
Brillouin zone was sampled at the T-point only and, un- 
less stated otherwise, the PBE XC functional has been 
employed 1 ^. Separable norm-conserving pseudopoten- 
tials were used to describe the interactions between the 
valence electrons and the ionic coresiii^ii 3 -. 

Long and well-equilibrated trajectories were necessary 
to obtain an accurate sampling. This requirement was 
made even more stringent by the strong dependence of 
the translational self-diffusion coefficient on temperature 
and, in the case of PBE water, on the expected low dif- 
fusivity at room temperatur o 103 ' 114 . Therefore, in each 
run, the system was equilibrated for 25 ps and the statis- 
tics were accumulated for the successive 250 ps. Finite- 
size effects are studied by comparing the results of the 
largest system with equally long runs on 64 and 32 water 



— 1 28 Water (PBE) 




FIG. 7. Comparison of goo(r) as obtained from second- 
generation CPMD simulations with 32, 64 and 128 water 
molecules. 

molecule. Two shorter 25 ps BOMD reference calcula- 
tions with 128 molecules were carried out using either 
Newtonian or Langevin dynamics to assess the accuracy 
of the simulations. The settings for both runs were iden- 
tical and started from the same well-equilibrated config- 
uration. The influence of the XC functional on the prop- 
erties calculated was investigated in a series of additional 
runs using a variety of different semilocal XC function- 
^86,90.115- 117, -p ne gt a tistics in each of these reference 
runs were accumulated for at least 30 ps after an equi- 
libration of 20 ps, resulting in more than 1 ns of AIMD 
simulations. 

All the simulations were performed using the 
CP2K/Quickstep code^ 3 -^, which relies on the mixed 
Gaussian and plane wave representation of the electronic 
degrees of freedom^. In this approach, the Kohn-Sham 
(KS) orbitals are expanded in terms of Gaussian orbitals, 
while a plane wave basis is used for the electron den- 
sity. Such a dual basis approach combined with advanced 
multigrid, sparse matrix and screening techniques en- 
ables one to achieve an efficient linear-scaling evaluation 
of the KS matrix. Efforts towards a full inear scaling 
algorithm are under wayi^-^ -. Here, the orbitals were 
represented by an accurate triple-^ basis set with two 
set of polarization functions (TZV2P) 5 -, while a density 
cutoff of 320 Ry was used for the charge density. 

Structural properties. The influence of finite-size 
effects on the oxygen-oxygen radial distribution function 
(RDF) goo( r ) is shown in Figure [7l It can be seen that 
the RDF is quickly converging with respect to system 
size: the goo( r ) for the 64-molecule system coincides 
within statistical uncertainties with the result of a larger 
128-molecule simulation. 

The accuracy of the second-generation CPMD can be 
readily established by comparing the calculated RDFs to 
the BOMD reference (Figure[8]). While the agreement be- 
tween these two methods is excellent, they both predict 
water to be overstructured, which can be seen by com- 
paring the present RDFs with the ones obtained from re- 
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cent neutron diffraction-^ and X-ray scattering^ exper- 
iments (Figure [5J . The most pronounced disagreement 
between the calculated and experimentally derived RDFs 
is the case of gonir), for which the relative heights of 
the first two intramolecular peaks is reversed. However, 
the inclusion of nuclear quantum effects 1 ^ 3 -*^, as well as 
London dispersion forcesi^Zr-i 3 ^ is expected to improve 
agreement with experiment. Using either artificially in- 
creased temperatures 1 ^ 2 - -10 ^ 132 or different XC function- 
a lg99-ioi ma y a i so j eac j ^ a better agreement with the 

experimental data. 

After verifying that the RDFs are converged with re- 
spect to the system size (i.e. finite-size effects are neg- 
ligible) a series of simulations with commonly used XC 
functionals were performed. Figure |H] and Table ITVl show 
that the shapes of the RDFs depend strongly on the par- 
ticular XC functional, whereas the coordination number 
for a water molecule as obtained by integrating the RDF 
curve up to the first minimum remains very close to four 
for all XC functionals. 

Obtaining a meaningful coordination number has im- 
portant implications for the debate around the structure 
and symmetry of the HB network in liquid water. The de- 
bate stems from an interpretation of the X-ray spectra of 
water as evidence for a large fraction of molecules (~80%) 
with broken HBs and a two-fold coordination. To check 
the validity of this interpretation, their spcctroscopically 
derived HB definition 2 - was used to calculate the average 
number of bonds formed by a water molecule. The re- 
sults presented in Table [V] show that the average coordi- 
nation number depends somewhat on the XC functional 
and varies between 3.1 (OLYP) to 3.6 (PBE). Given the 
fact that these results do not change qualitatively upon 
altering the HB dcfinition 128 i 134 i 135 , our AIMD simula- 
tions support a conventional nearly-tctrahedral view on 
the structure of liquid water. In the next section, we will 
analyze the symmetry and average number of HBs us- 
ing their electronic signatures in addition to the simple 
geometric criteria utilised here. 

Dynamical properties. Compared to the structural 
properties, the translational diffusion constant is known 
to exhibit stronger dependence on the size of the simu- 
lation box. The former is due to the fact that a diffus- 
ing particle entails a hydrodynamic flow, which decays 



rather slowly as 1/r. In a periodically repeated system 
this leads to an interference between a single particle and 
its periodic images. Analysing this effect, Diinweg and 
Kremer-i 3 ^ have derived the following relation: 



D(oo) = D(L) 



k B T( 
6ttt]L ' 



(19) 



where D(L) is the translational diffusion coefficient cal- 
culated for a system with the length of the simulation 
cell L, rj is the translational shear viscosity, whereas the 
constant £ is equal to 2.837. 

Figure [TO] shows that the 1 /L dependence predicted 
by Eq. [TO] holds rather well for the present AIMD sim- 
ulations, and thus, can be used to determine extrap- 
olated values for the translational diffusion coefficient 
D(oo) = 0.789 x 10~ 5 cm 2 /s and the shear viscosity 
coefficient r\ = 21.22 x 10~ 4 Pa • s. Comparison with 
the experimental values D = 2.395 x 10 _5 cm 2 /si2I and 
rj = 8.92 x 10~ 4 Pa-9i2£ confirms that liquid water as ob- 
tained from PBE AIMD simulations is less fluid than real 
water. However, taking into account that nuclear quan- 
tum effects were completely neglected in the simulations, 
it is possible to conclude that the PBE XC functional pro- 
vides a rather good representation of liquid water, much 
better than generally appreciated. 

Furthermore, Figure [TO] shows that applying Eq. [TO] 
together with the calculated r\ leads to an estimate of 
D(po), which is consistent with the extrapolated value. 
This demonstrates that 77 is much less system size de- 
pendent than D(po), and that the Stokes-Einstein rela- 
tion, which predicts an inverse relation between these two 
quantities, is no longer valid on nanometre scale. 

The velocity-velocity autocorrelation function for the 
hydrogen and oxygen atoms, as well as its temporal 
Fourier transform that represents the vibrational density 
of states, is shown in Figure[TT] The latter is of particular 
interest, since, beside being in excellent agreement with 
our BOMD reference calculations, it provides informa- 
tion about the dynamics of the HB network in liquid wa- 
ter. The small shoulder of the peak in the high-frequency 
oxygen-hydrogen stretching band indicates only an in- 
significant presence of dangling HBs, which implies that 
the time-averaged charge distribution in liquid water is 



TABLE IV. The position and height of the first maximum 
and minimum of goo{r) and the coordination number N c as 
calculated by integrating goo{r) up to the first minimum. 



XC 


max / \ 

9oo M 


max 

r oo 


ffooM 




N c 




3.25 


2.73 


0.44 


3.28 


4.04 


RPBB±i£ 


3.19 


2.75 


0.42 


3.32 


4.03 


revPBB±i£ 


3.01 


2.77 


0.50 


3.31 


4.05 


BLYP2&22 


2.92 


2.79 


0.57 


3.33 


4.09 


OLYP 2°iiI 


2.57 


2.79 


0.71 


3.30 


3.90 


Sopei^i 


2.75 


2.73 


0.78 


3.36 




ALSi22 


2.83 


2.73 


0.80 


3.4 


4.7 



TABLE V. The relative occurrence of double donor (DD), 
single donor (SD), no donor (ND) water molecules, percentage 
of donated and free HBs, as well as the average number of HBs 
formed by a water molecule. The results were obtained from 
AIMD simulations using several semi-local XC functionals. 



PBE 


RPBE revPBE 


BLYP 


OLYP 


DD 82.8 % 
SD 16.6 % 
ND 0.7 % 


81.4 % 
17.8 % 
0.8 % 


76.8 % 
22.0 % 
1.2 % 


72.9 % 
25.4 % 
1.7 % 


59.0 % 
36.3 % 
4.7 % 


donated HB's 91.0 % 
free HB's 9.0 % 


90.3 % 
9.7 % 


87.8 % 
12.2 % 


85.6 % 
14.4 % 


77.1 % 
22.9 % 



mean HB's 3.642 3.613 3.513 3.423 3.085 
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FIG. 8. Partial RDFs of liquid water at ambient conditions and its mean square displacement (top left panel). In the 
corresponding inset the onset of the cage effect can be observed at ~ 250 fs followed by diffusion, which is in excellent 
agreement with Gallo et alM^. 
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FIG. 9. Comparison of goo(r) obtained from neutron diffrac- 
tion, X-ray scattering and second-generation Car-Parrinello 
simulations using a variety of different XC functionals. 



mainly symmetric and the coordination is distorted but 
tctrahedral. 

Hydrogen bond kinetics. The kinetics of HB 
rearrangements is studied using the Luzar-Chandler 
modeli2^, which is able to describe the complex non- 
exponential relaxation behaviour of HBs with just two 
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FIG. 10. The translational diffusion coefficient as a function 
of system size computed by second-generation Car-Parrinello 
simulations using the PBE XC functional. The solid line is 
the linear extrapolation to D(oo) whereas the dotted line is 
the mean D(oo) obtained from Eg. 1191 



rate constants k and k' and the following reactive flux 
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FIG. 11. Velocity- velocity autocorrelation function and the 
corresponding vibrational density of states. The full lines are 
obtained from second-generation Car-Parrinello simulations, 
whereas the dashed lines represent BOMD reference calcula- 
tions. 



correlation function: 
dc(t) 



k(t) = 



((dh/dt) t=0 [1 - h(t)}) 



dt 
kc(t) - 



(h) 



k'n(t), 



(20) 



in which c(t) — (h(0)h(t)) / (h) is the HB autocorrelation 
function, n{t) = (h{0) [1 - h(t)} H{t))/(h), where H(t) is 
unity if the molecules of a selected pair are closer than 
R = 3.5 A and zero otherwise, while brakets (■) denote 
temporal averages. Thus n(t) represent the number of 
initially bonded pairs, that are broken at time t while 
remain closer than R. The HB population ope rator h(t) 
is defined using geometric descriptors from Ref. Il34 The 
HB lifetime is related to k by thb = k^ 1 , whereas the 
HB relaxation time is computed as 



/ dttc(t) 
Jdtc(t) ' 



(21) 



The reactive flux correlation function k(t) is non- 
cxponcntial and monotonically decaying after a period 
of few librations. A least squares fit of the simu- 
lation data obtained with the second-generation Car- 
Parrinello method to Eq. ((5DJ) yields fc=0.143 ps -1 and 
fc'=0.389 ps -1 , thus thb=6.98 ps and T r =10.25 ps. Com- 
pared to the results obtained with empirical poten- 
tials^, our AIMD values for thb and r r are both about 
twice as large, whereas the ratio 7v/thb=1.47 is very 
close to the value ~ 1.5 reported by other o 106 ' 139 . In addi- 
tion to these quantitative differences, c(t) obtained from 
ab initio simulations decays significantly slower than pre- 
dicted by he exponential law (Figure IT2"]) . In fact, we sus- 
pect that the decay might be biexponential, which would 
assume a second linear equation for n(t). Most likely this 
behaviour can be attributed to polarization as well as co- 
operativity effects, which are suposedly better described 
by DFT. 
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FIG. 12. The HB autocorrelation as a function of simulation 
time. 



Summary and further applications. The com- 
putationally efficient second-generation Car-Parrinello 
method extends AIMD simulations with hundreds of wa- 
ter molecules to the nanosecond timescale. These large- 
scale simulations allow the re-examination of the con- 
tribution of finite-size errors in the calculated structural 
and dynamical properties of liquid water, as well as the 
re-assessment of the performance of several commonly 
used XC functionals. It was found that structural prop- 
erties are well-converged and reproducible. By contrast, 
dynamical properties are much less established and the 
novel simulation approach enables the first calculation of 
the shear viscosity and the HB relaxation time of liquid 
water from first principles. It is demonstrated that AIMD 
simulations with classical nuclei and the PBE XC func- 
tional predict water to be overstructured and less fluid 
than the real water, with a shear viscosity within a factor 
of three from the experimental value. Given the fact that 
nuclear quantum effects are not included in the simula- 
tions it can be concluded that PBE provides a qualitative 
realistic model for interactions between water molecules 
in the liquid phase. 

In addition to the study reviewed here, the second- 
generation CPMD has been of great use in a recent in- 
vestigation of the structure of the water-air interface^ 
and for the calculation of absolute thermodynamic func- 
tions of liquid water from first principles'^. 



B. Electronic signature of the instantaneous 
asymmetry in the first coordination shell of liquid 
water 

As just presented, our AIMD simulations are consis- 
tent with the conventional wisdom that the local struc- 
ture of liquid water at ambient conditions is tetrahc- 
dralt^ 3 -. Although the thermal motion causes distor- 
tions from the perfectly tetrahcdral configuration, each 
molecule in the liquid is bonded, on average, to the 
four nearest neighbours via two donor and two acceptor 
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bondsi 4 ^. This traditional view is based on results from 
X-ray and neutron diffraction experiment s 121 ' 122 ' 144 , vi- 
brational spectroscopy^^ 8 -, macroscopic thermody- 
namics dat a 142 ' 143 ' 149 as well as molecular dynamics sim- 
ulations 1 ^^^^^^. 

Nevertheless, as mentioned above, this traditional pic- 
ture has recently been questioned based on data from the 
X-ray absorption, X-ray emission and X-ray Raman scat- 
tering experiments^iSr-i^ 4 -. The results of these spec- 
troscopic studies have been interpreted as evidence for 
strong distortions in the HB network with highly asym- 
metric distribution of water molecules around a central 
molecule. It has been suggested that a large fraction of 
molecules form only two strong HBs: one acceptor and 
one donor bond 2 -*!^ 2 --^. However, the "rings and chains" 
structure of liquid water—, as well as the inhomogeneous 
two-state mode l 152 ' 157 implied by such an interpretation, 
have been challenged on many front s 143 ' 147 ' 158 ' -^ and 
are a matter of an ongoing debate^ - ' 143 ' 147 ' 155 ' 164 ' 166 ' 172 . 

In this section, we review a computational study of 
the energetics and symmetry of local interactions be- 
tween water molecules in the liquid phase performed with 
ALMO EDAii 3 -. As demonstrated for the water dimcr, 
the decomposition of the interaction energy into physi- 
cally meaningful components provides a deeper insight 
into the nature and mechanisms of intermolecular bond- 
ing than the traditional total-energy electronic structure 
methods. Applying ALMO EDA to liquid water reveals a 
significant asymmetry in the strength of the local donor- 
acceptor contacts. DFT-based AIMD simulations per- 
formed using the second-generation Car-Parrinello ap- 
proac h 14 ' 108 enable us to characterise the geometric ori- 
gins of the asymmetry, its dynamical behaviour, as well 
as the mechanism of its relaxation. Furthermore, to ad- 
dress the controversial question of whether it is correct to 
interpret the X-ray spectra of water in terms of asymmet- 
ric structures, we present extensive calculations of its X- 
ray absorption (XA) spectrum and compare the spectral 
characteristics of water molecules with different degree of 
asymmetry. 

Computational details. ALMO EDA was per- 
formed for the dccorrelatcd configurations collected from 
a long 70 ps AIMD trajectory, which was generated by 
taking advantage of the computational efficiency of the 
second-generation Car-Parrinello method. AIMD sim- 
ulations with classical nuclei based on the PBE XC 
functionali^ - were performed at constant temperature 
(300 K) and density (0.9966 g/cm 3 ) for a system contain- 
ing 128 light water molecules. The results of the previous 
section show that PBE provides a realistic description of 
many important structural and dynamical characteristics 
of liquid water, including the RDFs, self-diffusion and 
viscosity coefficients, vibrational spectrum, and HB life- 
time. However, the PBE water with the classical descrip- 
tion of the nuclei is somewhat overstructured suggesting 
that the degree of the network distortion and fraction of 
broken bonds may be slightly underestimated. Therefore, 
we have explicitly verified that the main conclusions pre- 



sented here are also valid for snapshots generated with 
simulations with quantum hydrogen nuclei and based on 
a water model, which rectifies the main shortcomings of 
the PBE functional (see Supplementary Information in 
Ref.Qjl. 

The MOs in our ALMO EDA calculations were rep- 
resented by a triple-C Gaussian basis set with two sets 
of polarization functions (TZV2P) optimized specifically 
for molecular systems^. A very high density cutoff of 
1000 Ry was used to describe the electron density. The 
XC energy was approximated with the BLYP XC func- 
tional 8 ^ 9 ^. The Brillouin zone was sampled at the T- 
point only and separable norm-conserving pseudopoten- 
tials were used to describe the interactions between the 
valence electrons and the ionic cores 1 -^-^. Performing 
ALMO EDA with the HSE06 screened hybrid XC func- 
tionali 7 - 4 -, which provides better description of band gaps 
and electron derealization effects than BLYRiZ^-i 7 - 8 -, 
does not change the main conclusions presented below 
(see Supplementary Information in ReL^ 7 - 3 -). 

The XA calculations were performed at the oxygen Pl- 
edge using the half-core-hole transition potential formal- 
ismiSr-iSi within all-electron Gaussian augmented plane 
wave formalism to density functional theor y 182 ' 183 . The 
BLYP XC functional and large basis sets (6-311G** for 
hydrogen and cc-pVQZ for oxygen atoms) were used to 
provide an adequate representation of the unoccupied 
MOs in the vicinity of absorbing atoms 1 ^!. A density 
cutoff of 320 Ry was used to describe the soft part of 
the electron density. The onset energies of the absorp- 
tion spectra were aligned with the first ASCF excitation 
energy obtained for each oxygen atom from a separate 
calculation with the same setup. The spectral intensities 
were calculated as transition moment integrals in the ve- 
locity form. To mimic the experimental broadening, the 
discrete lines were convoluted with Gaussian functions of 
0.2 eV width at half- maximum. The final spectrum was 
obtained by averaging the convoluted spectra of 9,024 
oxygen atoms from 141 AIMD snapshots separated by 
500 fs (i.e. 64 oxygen atoms were randomly selected in 
each snapshot). 

Electronic asymmetry and its origins. Two- 
body derealization energy components AE x ^ y defined 
in Eq. are the main focus of this study. Each of these 
terms provide an accurate measure of the perturbation 
of orbitals localized on a molecule by donor or accep- 
tor interactions with a single neighbour in the bulk sys- 
tem. The donor-acceptor energies are calculated to in- 
clude cooperativity effects, which are essential for a cor- 
rect description of the HB networ k 40 ' 184 ' 185 . Further- 
more, the two-body terms are natural local descriptors 
of intermolecular bonding, which allowed us to analyse 
the molecular network in liquid water without introduc- 
ing any arbitrary definitions of a HB 1 ^. 

The electron derealization energy per molecule AEq 
can be analysed by neglecting the small higher order re- 
laxation term AEho (Figure fTB")) and by considering each 
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FIG. 13. Average contributions of the five strongest acceptor 
(AEn^c) and donor (AEc^n) interactions. The rightmost 
column shows that higher-order derealization (AEho) does 
not significantly contribute to the overall binding energy. 



water molecule either as a donor or as an acceptor: 



Mol 



Mol 



AEr 



&ec^n = AE 



(22) 



JV=1 



where C stands for the central molecule and N for its 
neighbours. It is important to emphasise that the terms 
donor and acceptor are used here to describe the role of 
a molecule in the transfer of the electron density. This is 
opposite to the labelling used for a donor and an acceptor 
of hydrogen in a HB. 

Figure [13] shows contributions of the five strongest 
donor-acceptor interactions to the average derealization 
energy of a molecule (AEc)- Brackets (. . .) denote av- 
eraging over all central molecules and AIMD snapshots, 
which were obtained by performing ALMO EDA for 701 
AIMD snapshots separated by 100 fs each (i.e. 89,728 
molecular configurations) . Figure H3] demonstrates that 
electron derealization is dominated by two strong inter- 
actions, which together are responsible for ^93% of the 
derealization energy of a single molecule. The third and 
the fourth strongest donor (acceptor) interactions con- 
tribute ~5% and correspond to back-donation of elec- 
trons to (from) the remaining two first-shell neighbours 
(i.e. there is non-negligible derealization from a typical 
acceptor to a typical donor). The remaining ~2% corre- 
spond to the derealization energy to (from) the second 
and more distant coordination shells. 

Comparison of the strengths of the first and second 
strongest donor-acceptor interactions (^25 kJ/mol and 
~12 kJ/mol, respectively) with that in the water dimcr 
(~9 kJ/mol) suggests that each water molecule can be 
considered to form, on average, two donor and two ac- 
ceptor bonds. Substantial difference in the strengths of 
the first and second strongest interactions implies that a 
large fraction of water molecules experience a significant 
asymmetry in their local environment. The asymmetry 
of the two strongest donor contacts of a molecule can be 



FIG. 14. The normalised probability density function of the 
asymmetry parameters Ta and Yd. The probability of find- 
ing a molecule in a bin can be found by dividing the corre- 
sponding density value by the number of bins (i.e. 400). The 
dashed black lines at T w |, |, | partition all molecules into 
four groups of equal sizes. 



characterised by a dimcnsionless asymmetry parameter 



T, 



1 



AE. 



AE, 



(23) 



An equivalent parameter Ta can be used for the two 
strongest acceptor contacts. The asymmetry parame- 
ter is zero if the two contacts are equally strong and 
equals to one if the second contact does not exist. The 
probability distribution of molecules according to their 
T-parameters is shown in Figure [14] together with the 
lines separating the molecules into four groups of equal 
sizes with different asymmetry. The shape of the dis- 
tribution demonstrates that most molecules form highly 
asymmetric donor or acceptor contacts at any instance 
of time. The line at T w 0.5, for example, indicates that 
for ^75% of molecules either or T r> is more than 
0.5, which means that the strongest donor or acceptor 
contact is at least two times stronger than the second 
strongest for these molecules. Furthermore, the peak in 
the region of high T in Figure [TJ] indicates the presence 
of molecules with significantly distorted or even broken 
HBs. 

Comparison of the configurations of donor-acceptor 
pairs involved in the first and second strongest inter- 
actions reveals the geometric origins of the asymme- 
try. It has been found that the strength of the inter- 
action is greatly affected by the intermolecular distance 
R = d(0 D - O a ) and the HB angle /? = ZO d O a H, 
while the other geometric parameters have only a minor 
influence on AEjj^a- The strongly overlapping distri- 
butions in Figure [15] suggest that some second strongest 
interactions have the same energetic and geometric char- 
acteristics as the strongest contacts. This implies that 
the electronic asymmetry observed cannot be attributed 
to the presence of two distinct types of HBs - weak and 
strong. It is, rather, a result of continuous deformations 
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FIG. 15. Distribution of the strength (A), intermolecular dis- 
tance 7? (B) and HB angle ft (C) for the first (red) and second 
(blue) strongest donor interactions C —> N. Filled areas show 
the contribution of configurations, for which back-donation to 
a nearby donor is stronger than donation to the second accep- 
tor (cf. text). Distributions for acceptor interactions, N — i> C, 
are almost identical and not shown. 



of a typical bond. Another important conclusion that 
can be made from the distributions in Figure [15] is that 
relatively small variations of the intermolecular distance 
(R ~ 0.2 A) and HB angle (/3 ~ 5 - 10°) entails remark- 
able changes in the strength and electronic structure of 
HBs. Analysis of the structure of the molecular chains de- 
fined by the first strongest bonds (i.e. one donor and one 
acceptor for each molecule) shows that their directions 
are random, without any long-range order (i.e. rings, 
spirals or zigzags) on the length scale of the simulation 
box (~ 15 A). 

The slow decay of the distribution tails (Figure [PS")) 
implies that it is difficult to quantify the concentration 
of single-donor and single-acceptor molecules in the liq- 
uid because defining such configurations using a distance, 
angle or energy cutoff is an unavoidably arbitrary proce- 
dure. A quantitative analysis of the HB network, which 
was performed in the previous section, shows that, ac- 
cording to the commonly used geometric definitions of 
HBs£ ' 134 i 135 , the structure of water is distorted tctrahe- 



drally with only a small fraction of broken bonds. How- 
ever, the results presented in Figure [T51 indicate that ge- 
ometric criteria cannot fully characterise the dramatic 
effect of distortions on the local electronic structure and 
donor-acceptor interactions of water molecules. 

It is important to note that some second strongest in- 
teractions are weakened by distortions to such an extent 
that back-donation to (from) a nearby donor (acceptor) 
becomes the second strongest interaction. Such config- 
urations can be clearly distinguished by the large-angle 
peak in Figure IT5TJ (the region of the 10-fold magnifica- 
tion). They account for ~ 6 — 7% of all configurations 
and are responsible for the low-energy peak in the distri- 
bution of AEd^-a (filled blue areas in Figure ITS"]) and for 
the high-T peak in Figure [HI 

Relaxation of the instantaneous asymmetry. 
The overlapping distributions in Figure [T3] suggest 
that, despite the difference in the strength of the 
donor-acceptor contacts, their nature is similar and 
the strongest interacting pair can become the second 
strongest in the process of thermal motion and vice versa. 
To estimate the time scale of this process, it is neces- 
sary to examine how the average energy of the first two 
strongest interactions fluctuates in time. The instanta- 
neous values at time t (solid lines in Figure ITBR) were 
calculated using the ALMO EDA terms for 3,501 snap- 
shots separated by 20 fs (448,128 local configurations) 
by averaging over time origins r separated by 100 fs and 
over all surviving triples: 



{AE c ^N(t)) 
1 T 



M(r,t) 

, MCT E AWt + *), (24) 



where M(r, t) is the number of triples that survived from 
time r to r+t. A triple is considered to survive a specified 
time interval if the central molecule has the same two 
strongest-interacting neighbours in all snapshots in this 
interval. 

Figure ITBR shows that the strength of the first two 
strongest interactions oscillates rapidly and after ~80 fs 
from an arbitrarily chosen time origin, the first strongest 
interaction becomes slightly weaker than the second 
strongest (note that first and second refer to their or- 
der at t = 0). The amplitude of the oscillations de- 
creases and the strengths of both interactions approach 
the average value of ^20 kJ / mol on the timescale of hun- 
dreds of femtoseconds. The decay of the oscillations in- 
dicates fast decorrelation of the time-separated instanta- 
neous values because of the strong coupling of a selected 
pair of molecules with its surroundings. In other words, 
although the energy of a particular HB fluctuates around 
its average value indefinitely (i.e. with a never-decreasing 
amplitude) , this bond has approximately equal chances of 
becoming weak or strong after a certain period of time in- 
dependently of its strength at t = 0. This effect is due to 
the noise introduced by the environment and can be ob- 
served in ultrafast infrared spectroscopy experiments^^. 
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The time averages shown in Figure [16] are physically 
meaningful and can be calculated accurately only for the 
time intervals that are shorter than the average lifetime 
of a HB thb ~ 5 — 7 ps, as shown in the previous sec- 
tio n 108 ' . The small residual asymmetry that is still 
present after 500 fs (Figure [TBR) is an indication of the 
slow non-exponential relaxation behaviour that charac- 
terises the kinetics of many processes in liquid water^i 

In addition to the instantaneous values of AEjj^a at 
time t, the dashed lines in Figure [ToT A show the corre- 
sponding averages over time t. These values were calcu- 
lated by averaging over time origins r, all snapshots lying 
in the time interval from r to r + 1 and over all surviving 
triples: 

(AE c ^ N {t)Y = 

1 T t M(r,t) 

= -V V— — - V AE c ^n(t + k)(25) 

Time-averages (A(t))* are related to the instantaneous 
values (A(t)) by the following equation: 

{A(t)r^- t j\A(r))dr, (26) 

where equality holds if all triples survive over time t. 

The dashed lines in Figure [ToT A shows that any 
neighbour-induced asymmetry in the electronic structure 
of a water molecule can be observed only with an exper- 
imental probe with a time-resolution of tens of femtosec- 
onds or less. On longer timescalcs, the asymmetry is 
destroyed by the thermal motion of molecules and only 
the average symmetric structures can be observed in ex- 
periments with low temporal resolution. 

An examination of the time dependence of all two-body 
and some three-body geometric parameters that charac- 
terise the relative motion of molecules reveals the mech- 
anism of the relaxation. Similar shapes of the curves in 
Figures \TE\k and [16)3 shows that the relaxation of the 
asymmetry is primarily caused by low-frequency vibra- 
tions of the molecules relative to each other. The minor 
differences in the behaviour of the curves, in particular 
at 80 fs, indicate that the relaxation of the asymmetry 
is also influenced by some other degrees of freedom. The 
temporal changes in the HB angles towards the average 
value fFigurcfToTJ) show that librations of molecules play 
this minor role in the relaxation process. 

The kinetics and mechanism of the asymmetry relax- 
ation presented here are supported by data from ultrafast 
infrared spectroscopy, which can directly observe inter- 
molecular oscillations with a period of 170 fsi^. They are 
also in agreement with the theoretical work of Fcrnandez- 
Serra et al., who have used the Mulliken bond order pa- 
rameter to characterise the connectivity and dynamical 
processes in the HB network of liquid watee^. 

X-ray absorption signatures of asymmetric 
structures. The time behaviour described above im- 
plies that the instantaneous asymmetry can, in princi- 
ple, be detected by X-ray spectroscopy, which has tem- 
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FIG. 16. Time dependence of the average strength (A), in- 
termolecular distance R (B) and HB angle /3 (C) for the first 
(red) and second (blue) strongest donor interactions C —¥ N. 
Solid lines show the instantaneous values while the dashed 
lines correspond to the time-average values. Time-dependent 
characteristics of acceptor interactions, N —¥ C, are almost 
identical and not shown. 



poral resolution of several femtoseconds and is highly 
sensitive to perturbations in the electronic structure of 
molecules2ii££. To identify possible relationships be- 
tween the spectroscopic features and asymmetry, the 
XA spectrum of liquid water is calculated at the oxy- 
gen K-edge. Although the employed computational ap- 
proach overestimates the intensities in the post-edge part 
of the spectrum and underestimates the pre-edge peak 
and overall spectral widths, it provides an accurate de- 
scription of the core-level excitation processes and semi- 
quantitatively reproduces the main features of the exper- 
imentally measured spectra (Figure [XTR ). 

The localized nature of the Is core orbitals allows the 
disentanglement of spectral contributions from molecules 
with different asymmetry. To this end, all molecules are 
separated into four groups according to the asymmetry of 
their donor and acceptor environments, as shown in Fig- 
ure [T3J Choosing boundaries at T « | , |, | distributes 
all molecules into four groups of approximately equal 
sizes (i.e. 25 ± 2%). Figure [TTB shows four XA spec- 



18 



tra obtained by averaging the individual contributions 
of molecules in each group. It reveals that molecules 
in the symmetric environment exhibit strong post-edge 
peaks, while molecules with a high asymmetry in their 
environment are characterised by the amplified absorp- 
tion in the pre-edge region. Furthermore, the relationship 
between the asymmetry and absorption intensity is non- 
uniform: the pre-edge peak is dramatically increased in 
the spectrum for the 25% of molecules in the most asym- 
metric group, for which the first strongest interaction is 
more than six times stronger than the second. As a con- 
sequence, the pre-edge feature of the calculated XA is 
dominated by the contribution of molecules in the highly 
asymmetric environments (Figure 1171 3). 

The pronounced pre-edge peak in the experimentally 
measured XA spectrum of liquid water has been inter- 
preted as evidence for its "rings and chains" structure, 
where ~80% of molecules have two broken HBg&l^. The 
results presented here suggest that this feature of the 
XA spectrum can be explained by the presence of a 
smaller fraction of water molecules with high instanta- 
neous asymmetry. Although the employed XA modelling 
methodology does not allow the precise estimation of the 
size of this fraction, this conclusion is consistent with 
that of recent theoretical studies at an even higher level 
of theory, which have demonstrated that the main fea- 
tures of the experimental XA spectra can be reproduced 
in simulations based on conventional nearly-tetrahedral 
modcl o 167 i 187 . Thus, the presented application of ALMO 
EDA to liquid water complements the previous results 
by revealing an interesting and important connection 
between relatively small geometric perturbations in the 
hydrogen-bond network, the large asymmetry in the elec- 
tronic ground state and the XA spectral signatures of the 
core-excitation processes. 

Summary. The main findings of our investigation of 
local donor-acceptor interactions of water molecules with 
their neighbours in the liquid phase are as follows: 

• The strength of donor-acceptor interactions suggest 
that each molecule in liquid water at ambient con- 
ditions forms, on average, two donor and two ac- 
ceptor bonds. 

• Even small thermal distortions in the tetrahedral 
HB network induce a significant asymmetry in the 
strength of the contacts causing one of the two 
donor (acceptor) interactions to become, at any 
instance of time, substantially stronger than the 
other. Thus, the instantaneous structure of water 
is strongly asymmetric only according to the elec- 
tronic criteria, not the geometric one. 

• Intcrmolccular vibrations and librations of OH 
groups of HBs result in the relaxation of the instan- 
taneous asymmetry on the timescale of hundreds of 
femtoseconds. 

• The pronounced pre-edge peak observed in the 
XA spectra of liquid water can be attributed to 
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FIG. 17. A. Calculated and experimental XA spectra of 
liquid water. B. Calculated XA spectra of the four groups 
of molecules separated according to the asymmetry of their 
donor (To) and acceptor (Ta) environments as shown in the 
inset. C. Contributions of the four groups into the total XA 
spectrum. The colour coding is shown in the inset above. 



molecules in asymmetric environments created by 
instantaneous distortions in the fluctuating but on 
average symmetric HB network. 



V. CONCLUSION 

We conclude by noting that the development of second- 
generation CPMDii as well as ALMO EDAi^^ are 
both for themselves important methodological achieve- 
ments that had been instrumental to elucidate the water 
dimer and liquid water, respectively. Nevertheless, nei- 
ther method alone but only the combination allowed for 
the additional insights to eventually reconcile the two ex- 
isting and seemingly opposite models of liquid water - the 
traditional symmetric and the recently proposed asym- 
metric - and represent an important step towards a bet- 
ter understanding of the electronic structure of the HB 
network of one of the most important liquids on Earth. 
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